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Abstract 

USM3D is a widely-used unstructured flow solver for simulating inviscid and viscous flows 
over complex geometries. The current version (version 5.0) of USM3D, however, does not 
have advanced turbulence models to accurately simulate complicated flow. We have imple- 
mented two modified versions of the original Jones and Launder k-e two-equation turbulence 
model and the Girimaji algebraic Reynolds stress model in USM3D. Tests have been con- 
ducted for three flat plate boundary layer cases, a RAE2822 airfoil and an ONERA M6 wing. 
The results are compared with those from direct numerical simulation, empirical formulae, 
theoretical results and the existing Spalart-Allmaras one-equation model. 

1 Introduction 

The unstructured-grid methodology offers some significant advantages compared to the tra- 
ditional structured-grid method for simulating flows over complex aerodynamic shapes. This 
is mainly due to the promise that the construction of an unstructured grid around complex 
configurations, such as aircraft, requires much less time than a comparable block-structured 
grid. Furthermore, local refinement of unstructured grids can be carried out more easily to 
improve the accuracy of the simulation. While more work remains to be done to fully realize 
their potential, much progress has been made in modeling complicated flows on unstructured 
grids (see Mavriplis [1] for a review). 

One important phenomenon for complex viscous flows is turbulence, which is difficult to 
simulate due to the existence of a wide range of scales. There are many types of methods 
to deal with turbulence, ranging from the simplest algebraic model to the most accurate 
direct numerical simulation. For most of the turbulence models, the Reynolds stress is as- 
sumed to be related to the mean strain rate by the eddy viscosity. Such turbulence models 
may be divided into zero-equation models (i.e., algebraic model), one-equation models, and 
two-equation models depending on the number of transport equations needed to be solved 
to obtain the eddy viscosity. The Reynolds stress model does not use the concept of eddy 
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viscosity. Instead, a transport equation for each component of the Reynolds stress tensor 
is solved directly. While the above models only solve the mean flow, large eddy simula- 
tions (LES) solve the large scale fluctuations in addition to the mean flow and only model 
the effects of small subgrid scales. Finally, the most accurate method is direct numerical 
simulation (DNS), where both the mean flow and all the fluctuations are solved directly. 

While the Reynolds stress model, LES, and DNS methods are much more accurate than 
the eddy viscosity based method, they require prohibitive amounts of CPU time and mem- 
ory. Therefore, the most widely used turbulence models in industry are still based on the 
concept of eddy viscosity. This is especially true for unstructured grid CFD codes. For ex- 
ample, the predominantly utilized turbulence model in the finite-element unstructured grid 
code of Mavriplis [2] is the algebraic model of Baldwin and Lomax [3], although extension 
has been made to include a two-equation model (Mavriplis and Martinelli [4]). The unstruc- 
tured nodal-based finite volume code FUN3D (see Anderson [5], Anderson and Bo n h a, ns [6]) 
contains two one-equation models, one by Baldwin and Barth [7] and the other by Spalart 
and Allmaras [8] . 

USM3D is a tetrahedral cell-centered unstructured flow solver for simulating inviscid 
and viscous flows over complex geometries. It was developed by Frink [9-12] at the NASA 
Langley Research Center and is now being widely used in both industry and government. 
This code is part of the TetrUSS flow analysis and design package which won the 1996 
NASA Software of the Year award. USM3D utilizes a cell-centered, upwind-biased, finite- 
volume, implicit /explicit algorithm capable of solving the compressible Euler and Navier- 
Stokes equations on unstructured tetrahedral meshes. Like most other unstructured CFD 
codes, however, USM3D does not have advanced turbulence models to accurately simulate 
complex flows. The current production version of the code has only the one-equation turbu- 
lence model developed by Spalart and Allmaras [8] , although a two equation k-e model was 
implemented in a previous research version of USM3D (Kwon and Hah 1 13]). It is well known 
that one-equation turbulence models are not adequate for complex flows such as separated or 
shear flows. The objective of this study is to present initial performance comparisons of the 
recently added advanced turbulence models with experimental data, theoretical results, DNS 
data, and the results from the existing Spalart-Alimaras one-equation model implemented 
by Frink [12]. 


2 Governing Equations 

The integral form of the Navier-Stokes equations, which govern the compressible Newtonian 
fluid flow in the absence of external forces, can be written as 

§,fff Q dV + JJ F(Q) • ndS = JJ G(Q )-ndS (1) 

where the first term (time changing rate) is integrated over the volume of a bounded do- 
main while the second term (convection) and the third term (diffusion) are integrated over 
the boundary of this domain. The quantities Q, F(Q) and G(Q) are vectors with five 
components and Eq. (1) contains five equations corresponding to the conservation of mass, 
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momentum, and energy. The unknowns in Eq. (1) are 
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where p, u, v. w and e are density, three velocity components, and total energy, respectively. 
The inviscid flux is 
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where p is the pressure and n x , n y , and n z are Cartesian components of the exterior surface 
unit normal n on the boundary of the domain. The viscous flux is 
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Moo is the free-stream Mach number and Rei is the Reynolds number based on a typical 
length scale (e.g., the total length of the flat plate). 

The total stress and heat flux g* can be divided into a laminar part (denoted by 
superscript L) and a turbulence part (denoted by superscript T) 


with the laminar part being 
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where p, L is the molecular viscosity, T is the temperature, Pr L is the molecular Prandtl 
number, and 7 = 1.4 is the ratio of specific heats. Following Eqs. (8) and (9), the Reynolds 
stress and turbulent heat flux can be approximated as 
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Note that the last term in Eq. (10) is presented only when a transport equation for turbulent 
kinetic energy k is solved. 


2.1 Transport Equations for k and e 

In the k-e model, two transport equations for the turbulent kinetic energy k and the dissi- 
pation rate e are solved. The eddy viscosity, / i T in Eqs. (10) and (11), is then calculated 
based on k and e 

P T = ( 12 ) 

where / M is a damping function and = 0.09. The transport equation for k is 
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Similarly, the transport equation for e can be written as 
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In all the equations as well as all the figures shown in this paper, unless stated explicitly 
otherwise, the length is normalized by a characteristic length L, the velocity by a oo, the 
density by p^, the viscosity by p^, the turbulent kinetic energy k by a^, and the dissipation 
rate e by p^ a^/poo, where a <*> is the free-stream speed of sound, p^ is the free-stream 
density, and p^ is the free-stream molecular viscosity. 
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The compressibility correction T in Eq. (14) and the damping function /), in Eq. (12) can 
take different forms. The two most widely used compressibility corrections are the Sarkar et 
al. [14] model 

T = M 2 (22) 

and the Wilcox [15] model 

r = (M 2 - Ml ') H(M t - M t o) (23) 

where H(x) is the Heaviside step function, the turbulent Mach number M t = \fkja with a 
being the local speed of sound. The damping function could take one of the following three 
forms: 

a) Jones and Launder [16] form 
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b) Van Driest form (Nagano and Hishida [17]) 


U = i - ex P (^r) ( 25 ) 

c) Speziaie et al. [18] form 

f “ = (‘ + w, ) tanh (^) (26) 

For the results shown in this paper, ivisc=6 refers to the Jones and Launder form of the 
damping function, Eq. (24), with no compressibility correction (i.e., F = 0) and ivisc=7 
refers to the Jones and Launder model modified by Carlson 1 19] as given in Eqs. (27-31). 
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The k-e model discussed above is the standard model, which is also called the “linear model” 
because the turbulent stress is linearly related to the mean strain rate by the eddy viscosity 
as is clear from Eq.(10). 
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2.2 Algebraic Reynolds Stress Model 

Various direct numerical simulation (DNS) data have shown that the turbulent stress does 
not vary linearly with the mean strain rate. A more accurate representation is the Reynolds 
stress model where a transport equation for each component of the Reynolds stress tensor 
is solved directly and the concept of eddy viscosity is not used. However, the Reynolds 
stress model requires a tremendous amount of CPU time and computer memory and thus, 
is seldom used in complex engineering applications. The algebraic stress model, in which a 
nonlinear term is added to the turbulent stress (so it is also called the “nonlinear model”), 
offers a practical compromise. Following Girimaji [20], the turbulent stress in the algebraic 
stress model is 
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where the mean strain rate and mean vorticity are 
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respectively. In Eq.(32), Ad and Ad are given by 
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C Q X = 3.4 C\ = 1.8 C 2 = 0.36 C 3 = 1.25 C A = 0.4 

(53) 

Furthermore, instead of = 0.09 for the linear model, = — G\ is applied 

model, which is invoked by setting ivisc=6 with inl=7 instead of inl=0. 

in the nonlinear 


3 Numerical Procedure 

The details of the numerical procedure for solving equation (1) are discussed in Frink [10-12] 
and only a brief overview is given here. The spatial computational domain is divided into 
a finite number of tetrahedral cells and a finite- volume discretization is applied to each cell. 
This procedure results in a set of volume-averaged state variables Q which are in balance 
with the area-averaged inviscid flux F(Q) and viscous flux G(Q). Inviscid fluxes are obtained 
across each cell face using either the Roe [21] flux- difference splitting approach or the Van 
Leer [22] flux-vector splitting technique. The data at nodal points could be obtained from the 
cell-averaged data by either an inverse-distance weighted averaging scheme or a Laplacian- 
weighted averaging scheme. The viscous fluxes are approximated at the cell-face centroids by 
linear reconstruction. An implicit time integration algorithm using the linearized backward 
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Euler time differencing approach is applied to update the solution. The resulting linear 
system of equations are solved at each time step with a subiterative procedure by a point- 
Jacobi method. Convergence to the steady state solution is accelerated by advancing the 
equations at each cell in time by the maximum permissible time step in that ceil. The CFL 
number is scaled according to the deviation of cell aspect ratio from the ideal value of an 
isotropic tetrahedron. 

The solution procedure for the two turbulence transport equations (13) and (17) are 
similar to that for solving the Navier-Stokes equations. Equations (13) and (17) are solved 
separately from the flow governing equations and from each other using the same backward 
Euler time integration scheme. The k and e equations can be solved using either first-order 
or second-order schemes. For the second-order method, either Roe’s SuperBee limiter or the 
Minimum Modulus (Min-Mod) limiter can be applied. To allow a large CFL number, an 
implicit time-stepping method is used to solve the k and e equations. 

Equations (13) and (17) are written in conservative form. There are two other forms 
of the transport equations: non-conservative form and incompressible form. For example, 
these three forms of the k equation can be written as 
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All three forms can be used to solve the k and e equations in USM3D (see § 4.1 for a 
comparison of the results). In general, the conservative form of Eq. (54) yields more accurate 
solutions than the other two forms, while the incompressible form of Eq. (55) is the least 
accurate since it is only applicable for incompressible flow. 

Two input parameters ko and p B are needed to specify the initial conditions and limit 
the smallest values of k and e. The initial conditions are k = ko and e = £o = C^pk^/ p B . 
The turbulent kinetic energy and dissipation rate are not allowed to become smaller than 
k Q and £o, respectively. On solid surfaces, the boundary conditions are k = k 0 and £ = 
L k / p{M 00 /Rei) 2 . Far field boundary conditions are applied by extrapolating k and £ from 
the interior for outflow boundaries and taken from the free-stream for the inflow boundaries. 

The k-e model has been coupled with wall function formulations to reduce the need for 
grid-resolving the flow in the near wall region. The following three different wall functions 
can be used together with the k-e model: (1) iwallf=l: the original wall function developed 
by Frink [12]; (2) iwallf=2: a wall function similar to that used in PAB3D (see Abdol- 
Hamid et al. [23] for details): (3) iwallf=3: same as the original wall function except the 
velocity components are zero at the solid surface. For iwallf=l, a slip velocity boundary 
condition is obtained by solving the Spalding formula 


y + = u + + e~ nB 


»+ _ i _ KU + _ ( KU+ ) 2 _ ( KU ' f 


(57) 
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with k = 0.4 and B = 5.5 using the Newton-Raphson iteration method while the velocity at 
the solid boundary for iwallf=2 and iwallf=3 is zero. 


4 Results 

4.1 Flat Plate Cases: BLT2, BLT3, BLT4 

Three tetrahedral grids for a simple flat plate boundary layer have been generated using the 
grid generator VGRID (Pirzadeh [24]). For each of the grids, results of the newly added 
k — g models are compared with the existing Spalart-Allmaras model and with experimental 
and theoretical data, (including DNS data) for the “wall function” grid, BLT2, and for the 
non-linear ARSM results on the near-wall resolved grids, BLT3 and BLT4. The computational 
domain extends from x = —0.5 to x 1.0 in the streamwise direction, from y = —0.02 to 
y = 0.02 in the spanwise direction for BLT3 and BLT4, (—0.05 < y < 0.05 for BLT2), and from 
z = 0 to z = 0.22 in the wall normal direction. The free-stream Mach number is = 0.5 
and Table 1 shows the Reynolds number for each case. In addition, Table 1 shows the grid 
size and minimum y + values. The grids for the flat plate cases are shown in Figures 1-3. 


Table 1: Summary of flat plate cases. 


Case 

Ceils 

Nodes 

Faces 

B. Nodes 

B. Faces 

Re L 

y + 

BLT2 

48,497 

9,629 

99,805 

2,813 

5,612 

2 x 10 e 

223 

BLT3 

37,483 

8,038 

78,328 

3,364 

6,724 

2 x 10 6 

1.8 

BLT4 

43,912 

9,321 

91,594 

3,772 

7,540 

20 x 10 6 

1.5 


There are four cells in the spanwise direction for BLT2 and two cells for the BLT3 and BLT4 
grids. Near the wall, the grid spacings for the BLT3 and BLT4 grids are much smaller than 
that for BLT2, but the grid spacings for BLT3 and BLT4 are much larger than that for BLT2 
far away from the wall. For BLT3 and BLT4, the first node point away from the wall has a 
y + & 1.8 and 1.5 respectively, while for the BLT2 grid y + « 223 at x/L = 0.5. Therefore, a 
wall function has to be applied for BLT2 while no wall function is needed for BLT3 or BLT4. 
Sample input files for the three plate cases are in Appendices A, B, and C. 

4.1.1 BLT2 Flat Plate - Wall Function Case 

The CFL number and residual history of the governing equations using the Spalart-Allmaras 
and linear k-e models are shown in Figure 4. The residual history for the Spalart-Allmaras 
and Carlson modified k — e model both converge rapidly at nearly the same rate to a level 
of 9 orders of magnitude smaller than the initial residual. The first linear model (ivisc=6) 
converged much slower after the 500 t/l time step and only decreased by less than 4 orders of 
magnitude after 4,000 time steps. The CFL number is allowed to increase dynamically from 
1 to 200 according to the residual: it increases when the residual is decreasing and decreases 
when the residual is increasing. For each of the three models tested, the CFL number rapidly 
increased to its maximum value within 40 time steps and remained at 200 for the duration 
of the run indicating a robust convergence. In Figure 5, velocity profiles are compared with 
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the empirical formula of Spalding, Eq. (57), while the skin friction coefficient is compared 
with the theoretical values for fully turbulent flow. 

C f = 0.0583 (Rea,) - 1/5 (58) 

It is evident that u + from the second (ivisc=7) k-e model is the closest to Spalding’s formula 
with the Spalart-Allmaras results the next closest. The first k-e model initially matches the 
data as dictated by the wall function, but then significantly under predicts u + . The skin 
friction coefficient from the first k-e model is significantly larger than the theoretical value 
whereas the Spalart-Allmaras model gives a result slightly smaller than the theoretical value 
and the second k-e model initially matches theory closely before predicting a slightly larger 
value. 

In summary, for the coarse flat plate grid using a wall function, (BLT2), the Carlson 
modified linear k-e model performed the best, with the Spalart-Allmaras model nearly as 
good and the first k-e model not in good agreement with empirical or theoretical data. 
The CPU time per timestep per ceil on an Intel Pentium II 300 MHz was 229/rsec for the 
Spalart-Allmaras model and 13% more for the linear k-e models. 

4.1.2 BLT3 Flat Plate - Grid Resolved Case 

The CFL number and the residual of the governing equations using the Spalart-Allmaras and 
linear k-e models are shown in Figure 6. The residual history for all three models converge 
rapidly at roughly the same rate to a level of 5 orders of magnitude smaller than the initial 
residual. The temporary hump in the residual trace around n = 1000 is characteristic of 
the establishment of turbulence for a non-wall function grid. For each of the three models 
tested, the CFL number rapidly increased to its maximum value within 75 time steps and 
remained at 200 for the duration of the run indicating a robust convergence. 

Inspection of the velocity profile and skin friction plots, see Figure 7, confirms that for 
the BLT3 grid a wall function is not needed. It is evident that u + from the Spalart-Allmaras 
model is closer to the Spalding curve with the second k-e model under predicting u + slightly 
and the first k-e model under predicting significantly. The skin friction coefficient from the 
first k-e model is significantly larger than the theoretical value whereas the Spalart-Allmaras 
model gives a result slightly smaller than the theoretical value and the second k-e model 
matches theory the closest predicting a slightly larger value. 

The turbulent kinetic energy and dissipation rate from the first k-e model (ivisc=6) at 
different downstream locations are compared with those from the DNS work of Spalart [25] 
in Figure 8. The turbulent kinetic energy k has been normalized by ( u *) 2 , where u* is the 
friction velocity, while dissipation rate e has been normalized by poo(u*) A / /i^. The profiles 
of both k and e at different downstream locations are very similar. The turbulent kinetic 
energy has a peak at y + ee 20 since k = 0 both at the wall and at the free stream. Although 
the peak value from the first k-e model is smaller than that from DNS, the overall agreement 
is good. There is no peak for the dissipation rate since e is the largest at the wall. The 
dissipation rate from the first k-e model also agrees well with that from DNS. 

The Reynolds stress components from the first k-e model are compared with those from 
the DNS of Spalart [25] in Figure 9. Similar to the turbulent kinetic energy, the Reynolds 
stress components have been normalized by (u*) 2 . The Reynolds shear stress from the first 


10 



k-e model agrees well with that from DNS. The three normal stress components from the 
first k-e model are almost identical whereas the streamwise component from DNS is much 
larger than the other two components. This is a property of the linear model, where the 
Reynolds stress is assumed to be isotropic. These three components will be different after we 
implement the nonlinear model (i.e., the algebraic stress model) by adding an extra no n li n ear 
term in the right hand side of equation (10). Note that the normal stress components, ww + , 
from k-e model in Figure 9(b) are very close to 2/3 k from Figure 8(a) because the first term 
in the right hand side of equation (10) is much smaller than the second term. Although not 
shown here, the three normal stress components from the Spalart-Allmaras model are nearly 
zero since the turbulent kinetic energy is not modeled (i.e., there is no 2/3 k term in equation 
(10) for this model). 

As discussed in § 3, the transport equations for k and e can be solved in three different 
forms: conservative, nonconservative and incompressible form. Figure 10 compares the 
residual history and streamwise velocity in wall coordinates from these forms. It is clear 
that the difference between the conservative and nonconservative forms is very small for this 
case. The residual from the incompressible form is much larger than the other two forms 
while the streamwise velocity from these three forms are similar. 

Although turbulent solutions can be obtained for the linear k-e models using an initial 
eddy viscosity of fi T / ' \i L = 1, it was found that this initial condition caused the boundary 
layer flow to remain laminar when used as the initial condition for the algebraic stress model. 
Therefore, the initial condition of /i T ///' = 10 was used for all algebraic stress cases. Figure 
11 shows the residual of the Navier-Stokes equations using the algebraic stress model with 
H 1 / n L = 10 for the flat plate boundary layer with the BLT3 grid. A great deal of oscillations 
occur and the residual settles down at about 0.75 orders of magnitude smaller than the initial 
residual after about 1,900 steps. Note that this residual is much larger than those from the 
first linear k-e model as well as the Spalart-Allmaras model in Figure 6. 

The skin friction coefficient and the streamwise velocity in wall coordinates from the 
nonlinear algebraic stress model are compared with those from the linear k-e model in Figure 
12(a). The streamwise velocity component from the nonlinear algebraic stress model agrees 
with Spalding’s empirical formula better than that from the first linear k-e model (ivisc=6). 
The skin friction coefficient from the algebraic stress model also has a better agreement with 
the theoretical value than that from the first linear k-e model from x ee 0.1 to x ee 0.5, at 
which point the linear model is more accurate. 

Since the success of the algebraic stress model depends on the initial condition, we ran a 
case where the algebraic stress model was restarted from the solution of the first linear k-e 
model. The skin friction coefficient and the streamwise velocity in wall coordinates from this 
run are shown in Figure 13. Similar to those shown in Figure 12, the streamwise velocity 
component agrees better with the Spalding empirical formula than that from the first linear 
k-e model while the skin friction coefficient also has a better agreement with the theoretical 
value than that from the first linear k-e model from x ee 0.1 to x ee 0.5, after which, the 
linear model is more accurate. 

In summary, for the BLT3 refined flat plate grid without using a wall function, the Carlson 
modified linear k-e model performed best for predicting skin friction and is a close second 
to Spalart-Allmaras for velocity profile. Experience with the linear models in the structured 
code PAB3D [23] indicates that with grid refinement the accuracy of the linear models 
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will improve significantly. The CPU time per timestep per cell on an SGI with two IP30 
processors was 194/rsec for the Spalart-Ailmaras model and 10% more for the linear k-e 
models. As discussed above, the residuals for the BLT2 and BLT3 grids using the algebraic 
stress model do not drop as much as those of the first linear k-e and the Spalart-Ailmaras 
models. Although both skin friction coefficient and streamwise velocity distributions from 
the algebraic stress model are better than those of the first linear k-e model, the results 
are still not as accurate as the Spalart-Ailmaras model. Furthermore, the results from the 
algebraic stress model depend on the initial condition. In an effort to decrease the streamwise 
distance needed for a turbulent boundary layer to develop, a high Reynolds number version 
of the BLT3 grid is examined in the following section. 

4.1.3 BLT4 Flat Plate 

The BLT4 grid was designed as a high Reynolds number counterpart to the BLT3 grid. It has 
a first nodal point normal distance of 1/10 of the BLT3 grid and was run with a Reynolds 
number is 20 million instead of 2 million. The CFL number and the residual of the governing 
equations using the Spalart-Ailmaras, both linear k-e models and the Girimaji algebraic 
Reynolds stress model are shown in Figure 14. The residual history for all cases except 
for the ARSM converge rapidly down 4 orders of magnitude, with the linear k-e models 
converging less robustly than the Spalart-Ailmaras model. In this case, the nonlinear model 
was run without using a linear solution as the initial condition, and despite setting the CFL 
number to a low constant value of 50 and running it for 8,000 iterations, the residual never 
settled down. The CFL number for the Spalart-Ailmaras and second linear models rapidly 
increased to its maximum value within 60 time steps while the first linear model initially 
increased then temporarily dipped to 25 around n = 120 before reaching the maximum 
allowable CFL at n = 160. 

Inspection of the velocity profile and skin friction plots, see Figure 15, confirms that for 
the BLT4 grid a wall function is not needed. It is evident that for this grid the Spalart- 
Ailmaras model is the closest to the experimental and theoretical data. The first linear 
model performs the worst and the second linear model and the ARSM are positioned between 
the other two, with the second linear model slightly closer to data than the ARSM. It is 
expected that with grid refinement and in less ideal configurations, the performance of the 
linear models and especially the nonlinear model will improve significantly. The CPU time 
per timestep per cell on an Intel Pentium II 300 MHz was 213/rsec for the Spalart-Ailmaras 
model, 13% more for the linear k-e models and 37% more for the non-linear model with 
nstagek=5. 

4.2 RAE2822 Airfoil 

In this section results for the previously mentioned turbulence models of Spalart-Ailmaras 
model, linear k—e (ivisc=6), and linear k — e (ivisc=7) will be compared with the Girimaji 
algebraic Reynolds stress model (ARSM), as well as, the experimental results of case 10 of 
Cook et al. [26] for a transonic airfoil. A sample input file is shown in Appendix D. 

The computational domain extends 6 chord lengths away from the airfoil in all directions 
with a channel width of 0.2c spanning 2 ceils. The free-stream Mach number of the boundary 
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layer, Reynolds number, and corrected free air angle-of-attack is M, x = 0.75, Re = 6.2 x 10 6 
and a = 2.81°, respectively. The grid contains 29,976 cells, 8,477 nodes, 66,772 faces, 6,820 
boundary nodes, and 13,640 boundary faces. The distance of the first node away from the 
wall corresponds to y + « 0.8. Full and closeup views of the grid are shown in Figure 16. 

Unlike the previous flat plate cases, the residual and CFL history (Figure 17) are much 
noisier, which may be due to unsteadiness in the region behind the shock. All of the models 
converge between 3 and 3.5 orders of magnitude after 4,000 time steps. Initially, the ARSM 
shows a great deal of oscillation since it is more dependent on the initial condition. Also, 
observed is a large increase in residual for the second linear model before finally settling 
down. The CFL history reveals that for the first 400 time steps each model is behaving 
similarly with the exception of the ARSM whose CFL history appears to be shifted to the 
left. After n = 400, the CFL number for all models oscillate periodically at high frequency. 

Lift and drag coefficient histories are shown in Figure 18. From the history plots, it 
is evident that the Spalart-Allmaras, first linear k-e, and Girimaji models follow the same 
convergence trend while the second linear model oscillates much more before settling down 
at around n = 2,500. The results at n = 4,000 along with experimental data and results 
from PAB3D for the first linear model are compared in Table 2. Final results show that 
the Spalart-Allmaras model is in closest agreement with experimental data. The next best 
results for USM3D are produced by the first linear model, the ARSM, and finally the second 
linear model. As an indicator of grid dependence for the advanced models, results in good 
agreement with experimental data from PAB3D are shown for the first linear model. 

In Figure 19, the predicted pressure coefficient from each model is plotted with the 
experimental data of Cook [26]. All models miss the location of the first suction peak by 
approximately 0.02c. This is likely due to the fact that the current implementation of the 
turbulence models in USM3D does not allow for setting a trip location for the boundary layer, 
while in the experiment it was set to 0.03c. Other areas of disagreement with experiment 
and among the models themselves is the shock location and region aft of the shock. The 
prediction of shock location by the second linear model was in excellent agreement while the 
other models all predicted different locations-all aft of the experimental location. Aft of the 
shock, all models predicted a lower Cp than experiment with the Girimaji model in the best 
agreement . 


Table 2: Comparison of Cl and Cp, with experimental data of Cook et al. [26] and PAB3D 
linear k-e model (ivisc=6) [23]. 


Case 

C L 

% diff 

c D 

% diff 

Experiment 

0.743 

- 

0.0242 

- 

Spalart-Allmaras 

0.732 

-1.5 

0.0214 

-11.6 

Linear k-e (ivisc=6) 

0.691 

-7.0 

0.0199 

-17.8 

Linear k-e (ivisc=7) 

0.649 

-12.7 

0.0178 

-26.4 

Girimaji ARSM 

0.803 

8.1 

0.0293 

21.1 

PAB3D Linear k-e (ivisc=6) 

0.720 

-3.1 

0.0257 

6.2 
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4.3 ONERA M6 Wing 

The tetrahedral viscous grid for the ONERA M6 wing was generated using VGRID and is 
similar in construction to those in Frink [12]. The grid contains 338,417 cells, 59,496 nodes, 
682,257 faces, 5,425 boundary nodes and 10,846 boundary faces. On the wing surface the 
midchord y + was approximately equal to 2 for the flow conditions of = 0.8447, Re mac » 
11.7 X 10 6 and a = 5.06°. The computational domain is bounded by a rectangular box defined 
by —6.5 < x < 6.5, 0 < y < 4, and —6.5 < 2: < 6.5, in aerodynamic coordinates relative to 
a semispan length of 1. Surface and symmetry plane meshes are shown in Figure 20. 

For this ONERA M6 wing case, results are presented for the Spalart-Allmaras model and 
the second linear (ivisc=7) model with varying numbers of Jacobi iterations in the solution 
of the k-e equations. The current implementation of the ARSM was found to be very stiff 
for the relatively coarse wing grid. Further investigations need to be performed to identify 
the minimum grid density required for a stable solution in various flow regimes. A sample 
input file is included in Appendix E. In Figure 21, the residual history is seen to be fairly 
similar among the models with the Spalart-Allmaras being the smoothest. Both linear cases 
are more noisy until n = 500. To aid solution stability, the original CFL ramping scheme in 
USM3D was used to bring the CFL number to a constant value of 75, rather than letting it 
vary dynamically up to 200 as in previous cases. 

In Figure 22, the coefficient of lift and drag histories confirm that the reduction in the 
number of iterations on the k-e equations causes no inaccuracy in the solution. This is 
significant because the reduction brings a 10% time savings. While the two linear cases 
agree with each other, the Spalart-Allmaras model predicts a 4% lower Cl while agreeing 
exactly with the linear models prediction of Co- 

Limiting surface streamlines simulating “oil-flow” patterns, shown in Figure 23, depict 
a significant amount of shock separated flow beyond the 77 et 0.65 for both models with the 
Spalart-Allmaras model results showing the strongest separation. The Cp plots in Figure 24 
demonstrate this with the exception of the 77 = 0.90 plane where Spalart-Allmaras matches 
the shock location well. It should be noted however, that the solutions in the tip region, 
77 > 0.90, can be particularly sensitive to a variety of factors such as grid density and 
turbulence model (see Rumsey and Vatsa [27]) and should be explored further in future 
studies. 

The CPU time per timestep per ceil on an SGI with two IP30 processors was 216 ysec 
for the Spalart-Allmaras model, 223 //sec for the second linear (ivisc=7) k-e model with 2 
Jacobi iterations (nstagek=2) on the k-e equations and 249 /isec for the second linear model 
with 6 Jacobi iterations. 

In summary, the linear k-e model is seen to perform well in a complex three-dimensional 
flow with only a small increase in CPU time required over the existing Spalart-Allmaras 
model. 


5 Concluding Remarks 

A systematic study has been conducted to assess the accuracy of two newly implemented 
turbulence models, modified versions of the standard linear two-equation k-e model and the 
non-linear algebraic Reynolds stress model of Girimaji. Initial results of test cases for the 
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flat plate, airfoil, and wing indicate that the new two-equation models yield comparable ac- 
curacy and efficiency to that of the Spalart-Ailmaras one-equation model. Work is currently 
underway to further examine factors such as grid sensitivities and solution details. 
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Figure 1: BLT2 grid for a flat plate boundary layer, where the near wall region is unresolved, 
a) whole domain: b) near wall region. 



Figure 2: BLT3 grid for a flat plate boundary layer, where the near wall region is resolved, 
a) whole domain: b) near wail region. 



Figure 3: BLT4 grid for a flat plate boundary layer, where the near wall region is resolved, 
a) whole domain: b) near wall region. 








Figure 4: a) Log residual: b) CFL history. BLT2 grid. 


u + vs. y + BLT2 Case Skin Friction BLT2 Case 




Figure 5: a) Streamwise velocity component in wall coordinates at x = 0.5: b) skin friction 
coefficient. BLT2 grid. 
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Figure 6: a) Log residual: b) CFL history. BLT3 grid. 


u + vs. y + BLT3 Case Skin Friction BLT3 Case 



Figure 7: a) Streamwise velocity component in wall coordinates at x = 0.5: b) skin friction 
coefficient. BLT3 grid. 
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a) 


b) 


Figure 8: Comparison of turbulent kinetic energy and dissipation rate from the first linear 
ks model (ivisc=6) with those from the DNS of Spalart [25]. a) turbulent kinetic energy; 
b) dissipation rate. BLT3 grid. 




a) b) 

Figure 9: Comparison of Reynolds stress components from the first linear k-e model 

(ivisc=6) at x = 0.5 with those from the DNS of Spalart [25]. a) shear component t ]2 ; 
b) normal components r^ a . BLT3 grid. 
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a) b) 

Figure 10: Comparison of solving the k and e equations in different forms using the first 
linear model (ivisc=6) at x = 0.5. (a) residual history; (b) velocity component in wall 
coordinates. BLT3 grid. 



Figure 11: Residual history of the governing equations using the algebraic stress model for 
the flat plate boundary layer solution on the BLT3 grid. 
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a) b) 

Figure 12: Comparison of a) streamwise velocity component in wall coordinate at x = 0.5 
and b) skin friction coefficient using the nonlinear algebraic stress and first linear k-e models 
for the flat plate boundary layer with the BLT3 grid. 




Figure 13: Skin friction coefficient and the streamwise velocity component at x = 0.5 in wall 
coordinates using the algebraic stress model (NL) restarted from the solution of the first 
linear k-e model (LN) and run for an additional 1100 and 2200 steps. BLT3 grid. 
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Log Residual for BLT4 Case 




b) 


Figure 14: a) Log residual: b) CFL history. BLT4 grid. 


u + vs. y + BLT4 Case Skin Friction BLT4 Case 



Figure 15: a) Streamwise velocity component in wall coordinates at x = 0.5 and b) skin 
friction coefficient for BLT4 grid. 
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Figure 16: RAE2822 airfoil grid, where the near wall region is resolved, a) whole domain; 
b) airfoil region. 


Log Residual for RAE2822 CFL History for RAE2822 
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n n 


a) b) 

Figure 17: a) Log residual: b) CFL history. RAE2822 airfoil. 
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C L History for RAE2822 



Experiment 
Spalart-Allmaras 
Linear k-e (ivisc=6) 
Linear k-e (ivisc=7) 
Girimaji ARSM 
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Spalart-Allmaras 
Linear k-e (ivisc=6) 
Linear k-e (ivisc=7) 
Girimaji ARSM 



4000 


a) 


b) 


Figure 18: a) Cl] b) Cd history. RAE2822 airfoil. 



Figure 19: Comparison of coefficient of pressure from each model at n = 4,000 with the 
experimental data of Cook [26]. 
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C L History OM6 Wing C D History OM6 Wing 



a) b) 


Figure 22: Coefficient of a) lift history: b) drag history. ONERA M6 wing case. 


OM6 USM3Dns Solution - 2000 cycles 

M=0.8447, AOA=5.06, Re=11.7M 



Figure 23: Comparison of k-e (ivisc=7) and Spalart-Allmaras model computed surface 
“oil-flow” patterns and flooded pressure contours for the ONERA M6 wing. 
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Cp 0M6 Wing r)=0.44 



C p 0M6 Wing ri=0.65 



C p 0M6 Wing r)=0.80 C p 0M6 Wing ri=0.90 



Figure 24: Comparison of coefficient of pressure from each model at n = 2, 000 with the 
experimental data of Schmitt [28]. 
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A BLT2 Input File 

USM3D (v5.0+turb) - Flat plate BL Second Linear Model with Wall Function 


Mach 

alpha 

beta 

ReUe ,mil 

Tinf ,dR 

itwall 

Tw/Tinf ipwall 

0.50000 

0.00000 

0.00000 

2.000000 

460.0 

0 

-1.0 0 

sref 

cref 

bref 

xmc 

ymc 

zmc 


1.00000 

1.00000 

1.0000 

0.00000 

0.00000 

0.00000 


impl 

dt/ cf 11 

iramp 

cf 12 

cf lmin 

nstage 

iautocf 1 

1 

-1.00 

4 

200. 

10.00 

5 

1 

irest 

mstage 

iresmth 

dqmax 

ptol 

frac_hit 

limiter 

0 

3 

1 

0.25 

.001 

1.0000 

0 . 

nupdate 

nwrest 

ipltqn 

if ast 

idiagnos 



10 

50 

1 

0 

1 



iorder 

lapl-avg 

high-bc 

ifds 

ivisc 

ckv 


2 

1 

1 

1 

7 

0.25 


ncyc 

nengines 

bcf ile 

compF&M 

cldes 

idim 


4000 

0 

1 

0 

0 . 

1 



Files: 1) Connectivity/grid, 2) Face/b.c., 3) Restart, 4) Flow output (a40) 
blt2.cogsg blt2.iface blt2.urest blt2.flo blt2.mapbc 
Begin KE Cont 

ilhg ikeord icons dtf itk icomp int ut/ul intb ut/ulb ini iwallf bke dkemax isk 
-14 1 0 1.0 2 0 l.e-2 10. l.e-2 10. 0 11. 0.00 0 

ratiokp nstagek 
0.00 5 

End KE Cont 

B BLT3 Input File 


USM3D (v5.0+turb) - Flat plate BL Second Linear Model without Wall Function 


Mach 

alpha 

beta 

ReUe ,mil 

Tinf ,dR 

itwall 

Tw/Tinf ipwall 

0.50000 

0.00000 

0.00000 

2.000000 

460.0 

0 

-1.0 0 

sref 

cref 

bref 

xmc 

ymc 

zmc 


1.00000 

1.00000 

1.0000 

0.00000 

0.00000 

0.00000 


impl 

dt/ cf 11 

iramp 

cf 12 

cf lmin 

nstage 

iautocf 1 

1 

-1.00 

4 

200. 

10.00 

5 

1 

irest 

mstage 

iresmth 

dqmax 

ptol frac_hit 

limiter 

0 

3 

1 

0.25 

.001 

1.0000 

0 . 

nupdate 

nwrest 

ipltqn 

if ast 

idiagnos 



10 

50 

1 

0 

1 



iorder 

lapl-avg 

high-bc 

ifds 

ivisc 

ckv 


2 

1 

1 

1 

7 

0.25 


ncyc 

nengines 

bcf ile 

compF&M 

cldes 

idim 


4000 

0 

1 

0 

0 . 

1 


Files: 1) Connectivity/grid, 2) 

Face/b . c . 

, 3) Restart 

4) Flow 

output (a40) 

blt3 . cogsg 

blt3 . if ace 

blt3.urest blt3.flo 

blt3 .mapbc 




Begin KE Cont 

ilhg ikeord icons dtf itk icomp int ut/ul intb ut/ulb ini iwallf bke dkemax isk 
-14 1 0 1.0 2 0 l.e-2 10. l.e-2 10. 0 01. 0.0 0 

ratiokp nstagek 
0.00 5 

End KE Cont 
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C BLT4 Input File 

USM3D (v5.0+turb) - Flat plate BL Second Linear Model without Wall Function 


Mach 

alpha 

beta 

ReUe ,mil 

Tinf ,dR 

itwall 

Tw/Tinf : 

0.50000 

0.00000 

0.00000 

20.00000 

460.0 

0 

-1.0 

sref 

cref 

bref 

xmc 

ymc 

zmc 


1.00000 

1.00000 

1.0000 

0.00000 

0.00000 

0.00000 


impl 

dt/ cf 11 

iramp 

cf 12 

cf lmin 

nstage 

iautocf 1 

1 

-1.00 

4 

200. 

10.00 

5 

1 

irest 

mstage 

iresmth 

dqmax 

ptol 

frac_hit 

limiter 

0 

3 

1 

0.25 

.001 

1.0000 

1 . 

nupdate 

nwrest 

ipltqn 

if ast 

idiagnos 



10 

50 

1 

0 

1 



iorder 

lapl-avg 

high-bc 

ifds 

ivisc 

ckv 


2 

1 

1 

1 

7 

0.25 


ncyc 

nengines 

bcf ile 

compF&M 

cldes 

idim 


4000 

0 

1 

0 

0 . 

1 



Files: 1) Connectivity/grid, 2) Face/b.c., 3) Restart, 4) Flow output (a40) 
blt4.cogsg blt4.iface blt4.urest blt4.flo blt4.mapbc 
Begin KE Cont 

ilhg ikeord icons dtf itk icomp int ut/ul intb ut/ulb ini iwallf bke dkemax isk 
-14 1 0 1.0 2 0 l.e-2 10. l.e-2 10. 0 0 1.0 0.00 0 

ratiokp nstagek 
0.00 5 

End KE Cont 

D RAE2822 Input File 

USM3D (v5.0+turb) - RAE2822 Airfoil Second Linear Model without Wall Function 


Mach 

alpha 

beta 

ReUe ,mil 

Tinf ,dR 

itwall 

Tw/Tinf ipwall 

0.75000 

2.81000 

0.00000 

6.200000 

460.0 

0 

o 

0 

H 

1 

sref 

cref 

bref 

xmc 

ymc 

zmc 


0.20000 

1.00000 

0.2000 

0.00000 

0.00000 

0.00000 


impl 

dt/ cf 11 

iramp 

cf 12 

cf lmin 

nstage 

iautocf 1 

1 

o 

0 

T—\ 

1 

4 

200. 

10.00 

5 

1 

irest 

mstage 

iresmth 

dqmax 

ptol 

f rac_hit 

limiter 

0 

3 

1 

0.25 

.001 

1.0000 

0 . 

nupdate 

nwrest 

ipltqn 

if ast 

idiagnos 



10 

50 

1 

0 

1 



iorder 

lapl-avg 

high-bc 

ifds 

ivisc 

ckv 


2 

1 

1 

1 

7 

0.25 


ncyc 

nengines 

bcf ile 

compF&M 

cldes 

idim 


4000 

0 

1 

0 

0. 

1 



Files: 1) Connectivity/grid, 2) Face/b.c., 3) Restart, 4) Flow output (a40) 
rae3.cogsg rae3.iface rae3.urest rae3.flo rae3.mapbc 
Begin KE Cont 

ilhg ikeord icons dtf itk icomp int ut/ul intb ut/ulb ini iwallf bke dkemax isk 
-14 1 0 1.0 2 0 l.e-2 1.0 l.e-2 1.0 0 0 1.0 0.25 0 

ratiokp nstagek 
0.05 5 

End KE Cont 
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E ONERA M6 Wing Input File 


USM3D (v5+turb) - ONERA M6 wing 


Mach 

alpha 

beta 

ReUe ,mil 

Tinf ,dR 

itwall 

Tw/Tinf ipwall 

0 . 84470 

5.06000 

0.00000 

21.191770 

460.0 

1 

-1.0 0 

sref 

cref 

bref 

xmc 

ymc 

zinc 


0.52550 

0.66700 

1.00000 

0.00000 

0.00000 

0.00000 


impl 

dt/ cf 11 

iramp 

cf 12 

cf lmin 

nstage 

autocf 1 

1 

-75.0000 

1 

75. 

10.0 

6 

0 

irest 

mstage 

iresmth 

dgmax 

ptol 

frac_hit 

limiter 

1 

3 

1 

0.25 

0.001 

0.001 

0 . 

nupdate 

nwrest 

ipltqn 

if ast 

idiagnos 

nodeypl 


10 

100 

2 

0 

2 

0 


iorder 

lapl-avg 

high-bc 

ifds 

ivisc 

ckv 


2 

1 

1 

1 

7 

0.25 


ncyc 

nengines 

bcf ile 

compF&M 

cldes 

idim 


500 

0 

1 

0 

0 . 

1 


Files: 1) Connectivity/grid, 2) 

Face/b . c . 

, 3) Restart, 4) Flow 

output (a40) 


om6ke_c . cogsg 
om6ke_c . if ace 
om6ke_c . urest 
om6ke_c . f lo 
om6ke_c.mapbc 
Begin KE Cont 

ilhg ikeord icons dtf itk icomp int ut/ul intb ut/ulb ini iwallf bke dkemax isk 
-14 1 0 1. 2 0 l.e-2 10. l.e-2 10. 0 01. 0.25 0 

ratiokp nstagek 
0.10 6 
End KE Cont 


31 



REPORT DOCUMENTATION PAGE 

Form Approved 
OMB No. 0704-0188 

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data 
sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other 
aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and 
Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), 
Washington, DC 20503. 

1 . AGENCY USE ONLY ( Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 

April 2000 Contractor Report 

4. TITLE AND SUBTITLE 

Implementation of Advanced Two Equation Turbulence Models in the 
USM3D Unstructured Flow Solver 

5. FUNDING NUMBERS 

WU 522-31-21-01 
NAS4-50066 

6. AUTHOR(S) 

Qunzhen Wang, Steven J. Massey, and Klialed S. Abdol-Hamid 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

Analytical Services & Materials, Inc. 

Lancaster, California 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23681-2199 

10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

NASA/CR-2000-210102 

11. SUPPLEMENTARY NOTES 

Langley Technical Monitor: Neal T. Frink 

12a. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified-Unlimited 

Subject Category 02 Distribution: Standard 

Availability: NASA CASI (301) 621-0390 

12b. DISTRIBUTION CODE 

13. ABSTRACT (Maximum 200 words) 

USM3D is a widely-used unstructured flow solver for simulating inviscid and viscous flows over complex 
geometries. The current version (version 5.0) of USM3D, however, does not have advanced turbulence models 
to accurately simulate complicated flow. We have implemented two modified versions of the original Jones and 
Launder k-epsilon "two-equation" turbulence model and the Girimaji algebraic Reynolds stress model in 
USM3D. Tests have been conducted for three flat plate boundary layer cases, a RAE2822 airfoil and an ONERA 
M6 wing. The results are compared with those from direct numerical simulation, empirical formulae, theoretical 
results, and the existing Spalart-Allnraras one-equation model. 

14. SUBJECT TERMS 

unstructured, tetrahedral, Navier-Stokes, turbulence model 

15. NUMBER OF PAGES 

36 

16. PRICE CODE 

A03 

17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 

18. SECURITY CLASSIFICATION 
OF THIS PAGE 

Unclassified 

19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 

20. LIMITATION 
OF ABSTRACT 

UL 


NSN 7540-01-280-5500 


Standard Form 298 (Rev. 2-89) 

Prescribed by ANSI Std. Z-39-18 
298-102 




























